{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c9e57fad-6e09-4f1b-b129-d1e500629f1a",
   "metadata": {},
   "source": [
    "# Structure relaxation\n",
    "\n",
    "This tutorial demonstrates how to relax a structure and calculate an energy-volume curve with only a few commands.\n",
    "\n",
    "First we need to import a few support functions for creating the structure (`bulk`), the `CPUNEP` calculator and the `relax_structure` function from calorine as well as functionality from `matplotlib` and `pandas` that we will use for plotting."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c783ad70",
   "metadata": {},
   "source": [
    "All models and structures required for running this tutorial notebook can be obtained from [Zenodo](https://zenodo.org/record/10658778).\n",
    "The files are also available in the `tutorials` folder in the [GitLab repository](https://gitlab.com/materials-modeling/calorine/-/tree/master/tutorials)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "96eba022-6823-44bf-856e-727e18d4e367",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ase.build import bulk\n",
    "from ase.units import GPa\n",
    "from calorine.calculators import CPUNEP\n",
    "from calorine.tools import relax_structure\n",
    "from matplotlib import pyplot as plt\n",
    "from pandas import DataFrame"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fca10a1-f6c9-4f08-b3c2-bfe5600ef0cb",
   "metadata": {},
   "source": [
    "Next we create a structure, set up a `CPUNEP` calculator instance, and attach the latter to the structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2bb0e155-fcf4-4edc-96e3-cb5b590c802d",
   "metadata": {},
   "outputs": [],
   "source": [
    "structure = bulk('PbTe', crystalstructure='rocksalt', a=6.6)\n",
    "calc = CPUNEP('nep-PbTe.txt')\n",
    "structure.calc = calc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8026639-040a-4e2b-9cb8-782843713689",
   "metadata": {},
   "source": [
    "Now we can relax the structure by calling the `relax_structure` function.\n",
    "To demonstrate the effect of the relaxation, we print the pressure before and after the relaxation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "069cfd9b-8792-4c03-893c-45f0d1df6038",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pressure before: 0.68 GPa\n",
      "pressure after:  0.00 GPa\n"
     ]
    }
   ],
   "source": [
    "pressure = -np.sum(structure.get_stress()[:3]) / 3 / GPa\n",
    "print(f'pressure before: {pressure:.2f} GPa')\n",
    "\n",
    "relax_structure(structure)\n",
    "pressure = -np.sum(structure.get_stress()[:3]) / 3 / GPa\n",
    "print(f'pressure after:  {pressure:.2f} GPa')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25cecc1c-a7f0-4ba9-8b1b-19911066d31b",
   "metadata": {},
   "source": [
    "We can readily calculate the energy-volume curve from 80% to 105% of the equilibrium volume.\n",
    "To this end, we use again the `relax_structure` function but set the `constant_volume` argument to `True`.\n",
    "This implies that the volume remains constant during the relaxation while the positions and the (volume constrained) cell metric are allowed to change.\n",
    "The results are first compiled into a list, which is then converted into a pandas `DataFrame` object.\n",
    "The latter is not necessary as we could simply use a list but it is often convenient to use `DataFrame` objects when dealing with more complex data sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e5376b3d-527a-4c67-b7a1-8a0c1852ad08",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = []\n",
    "for volsc in np.arange(0.8, 1.05, 0.01):\n",
    "    s = structure.copy()\n",
    "    s.cell *= volsc ** (1 / 3)  # the cubic root converts from volume strain to linear strain\n",
    "    s.calc = calc\n",
    "    relax_structure(s, constant_volume=True)\n",
    "    data.append(dict(volume=s.get_volume() / len(structure),\n",
    "                     energy=s.get_potential_energy() / len(structure),\n",
    "                     pressure=-np.sum(s.get_stress()[:3]) / 3 / GPa))\n",
    "df = DataFrame.from_dict(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35b7c27b-66ee-43fa-a8f3-e0e51af05c5a",
   "metadata": {},
   "source": [
    "Finally, we plot the energy-volume curve along with the pressure-volume curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "28199c21-2839-4a7c-91f3-fa70fd83adfc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEwCAYAAADB3LQ7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA3XAAAN1wFCKJt4AABkVklEQVR4nO3dd3xTVRvA8V/SRSfpgJbSwZQhUKbsvWSKAxBlD2W7FVEQFZRXZSiCoCggyhBZhSIKAjIsKJQCsmQVKKUto3u3ue8fsYFCW5KSNmn6fD+fKrn35uTp5XKenHvPUCmKoiCEEEJYKbW5AxBCCCGKkyQ6IYQQVk0SnRBCCKsmiU4IIYRVk0QnhBDCqkmiE0IIYdUk0QkhhLBqkuiEEEJYNVtzB1BaqFQqc4cghBDFzhrnEJFEZwRrvACEECKXtX6hl1uXQgghrJokOiGEEFZNEp0QQgijXbt2jX79+uHp6YmXlxcDBgzgxo0b5g4rX5LohBBCGG3ChAkAXL58mUuXLpGens7kyZPNHFX+JNEJIYQw2sWLFxkwYAAuLi64uroycOBATpw4Ye6w8iW9Lk1tz2y0LpUItunCytAIIuPT8NM4MqRlFfrm7ESdfB06TDF3lEII8VBeffVV1q1bR69evVAUhdWrV9OnTx9zh5UvadGZmNalEmx9mUPr5xF2JZ6YxAzCrsRzaP082Pqybr8QQlgolUql/5kxY0aBx7Vu3ZrY2Fjc3d3x8PAgLi6Ot99+u+QCNYK06Ews2KYLh7JHMdNmKc5KKqeUQPxVscy0+ZZ3skfT3KYL/cwdpBBCFMCQ8cJarZauXbsyYMAAduzYAcCMGTPo1q0bBw8eLO4QjaZSZBS0QVQqlUEXwNOLDhB2JZ4BNrv52PYb1CpQFFiZ04UZ2cNpFOjJ+nGtSiBiIYQwjqH13M2bN6lQoQJXr17Fz88PgKtXrxIQEMCNGzfw8vIq7lCNIrcuTSwyPg0FWJvTkQjFBwCVCoba7mS7/VvUv/kL5GSbN0ghhHgIXl5e1KhRg4ULF5Kenk56ejoLFy7Ez8/P4pIcSKIzOT+NIypgoM1uAlUxHNVWR6uoiNB684j6GjNyvoAFjeHvbyErHfbMhiMr8i/syArdfiGEsDCbN28mLCyMypUrU6lSJf766y+Cg4PNHVa+5NalgQxt0m86eo1D6+cx02YpU7NHszanIwNtdvOR7VKW5vSin2883rH7dQe7eIN/CzizBXrPhybD7hR0ZAVsffn+7UIIUUwMredKG0l0BjL0AtAe1iWod7JHsSa7I7nvyE122T3n4eDfEPbNhdNbAAVsHSE7Hbp/BC3HS5ITFkNRFLKzs62y8rN0KpUKW1vbEp1oWRJdGWfwBXD3OLqDl7kWl0pld0fSs7TUi9lM7yrQdsxnumNvnIX98+H4WlBydNtcvCHlhiQ5YXbZ2dlcunSJ7Gx5pmwuKpVKP8VWSSQ8SXRl3MNeAFHxaXSd+wfp2Vq2TGxDXV+3OzvjLsOfC+DvpYACqKDxEGj9MnhWf9jQhSiSyMhIMjMzqVy5stUu32LpUlNTiY6OxtPTkwoVKhT750miK+NMcQEsP3CJGVtOEeRXng3jW2OjvqvyyL1d6VIRkqL/+1A11O0HbV6BSg0e6rOFMIZWq+Xs2bMEBATg7Oxs7nDKtPj4eKKjo6lVq1axf+Gw1kQnvS5L0JCWVWjor+FYZALfh0bc2XH3M7nXzsLjswEV2DnByQ2wpC388Axc/lN6aYoSkZOju5VuZ2dn5kiEk5OT/lmpKBpJdCXIRq3i46fqY6tW8emvZ7kWn5Z/x5MW46DP55CVCvWegfL+cH4HLOuhe5639aX7k11uOa4yxZh4eNb4rb60k7+TopNEV8LqVHLjhXbVSM3MYdqmf1ASo/LveNJkmG67V02YfBT6fQVeteD2Rd1UK1smQ/Bk0OZIL00hhCiEPKMzkCnvXadn5dB9/l4u30pl4XON6dXAwFaYVgtnQ3RDE6LCdNtsHCAnE3rNgWajTBKfEJmZmVy4cIHq1atjb29v7nDKtJL8u5BndMJkytnZ8NGT9QF4L/gkCalZhr1RrYY6fWDMLhi6GRxcIScDUGDvZxC6CDJTii9wIQyg1SpsOnqNpxcdoPlHO3l60QE2Hb2GVmt9FagoHSTRmUnrGl4808SPm8kZzN5+2rg3q1S6IQmZKVChjm5bUhT8+jbMqwd/fAppcaYPWogH0GoVJq0O47V1x/IsU/XaumNMWh1WapNdVpaBX0aFRZJEZ0bv9KyDh7M9q/+6yqGLtwx/493P5CYchD5f6IYiVG4K6fGweybMqw873oPfpkkvTVFigo9Fsf1kDDlaRT8rkALkaBW2n4wh+FiUST8vOTmZiRMnEhAQQMWKFRk6dCgJCQlERESgUqlYuXIlNWrUQKPRMHz48DwJKywsjI4dO+Lh4UGNGjX45ptv9PtmzJhB7969GTduHB4eHkyZMoWMjAzGjh2Lh4cHVatW5dtvv0WlUhEREcGxY8dwdXUlOTlZX8a1a9dwcHAgKsq0v7MwnqxHZ0buzvZM712Xl9eGM3H1UfzcHYm6e0XyIF/U6nvGzeTX8ST3/1tfho7vQEIkhP8IB+aD2lbXYSXttm48Xn7lCGGA1rN3kZheeMsmJSObghptOVqFV38KZ9rmfwotw62cHQemdDIoppEjR2Jra8vx48exs7Nj9OjRTJw4kQ8//BCAX375haNHj5KUlETz5s358ccfGT58ONHR0XTt2pWvvvqKp59+mtOnT9OtWzeqVatG586dAdi+fTtLly5lwYIFZGZmMnPmTA4fPszJkydxcnLi+eef18cRFBRErVq1+Pnnnxk+fDgA33//PV26dMHX19eg30UUH2nRmVmfBpXwcnHgRlIGRw251ZN0vfBemtps6DMfXjoOLSfqOqugwM4ZsLSrbtox6aUpismD+jGYsp/DjRs3WL9+PQsXLkSj0eDs7MwHH3zA2rVr9eMAp0+fjqurK76+vjz++OMcOXIEgJUrV9KuXTsGDBiAjY0N9erVY8SIEaxatUpffr169Rg+fDi2trY4OTmxatUqpkyZQqVKlShfvjzvvfdennhGjRrF8uXL9a9XrFjBiBEjTPcLiyKzyBZdWFgYY8aM4dKlS2i1WurWrcvs2bNp165dge+JjIzklVde0a9227x5c3799VcAli9fzqhRo3B0dNQfP336dN58883i/UUMsOX4dW6nZOTZdu+tnn6NKt/Z2WFKwYXdnbTcKkH3WdD2NTi0BA58DpF/wcLHdPvbvyVJThjFkFZW7sLD+eUzFdA40N1kCw9HRESg1WqpWrVqnu1qtZroaN3sQj4+Pvrtzs7OxMfH69+7bds2NBqNfn9OTg5t27bVvw4ICMhTblRUFP7+/gXuHzRoEK+//jqXLl0iOjqamzdv0rdv34f6HYVpWGSiCwwMZMOGDfoLaePGjfTq1YvY2Ng8ySpXSkoKHTt2ZNiwYXz33Xc4Ojpy9OjRPMfUr1+f8PDwkgjfKCtDIwr8lqvVKqw8eDlvojOWkwd0fBtaTYRFLXS3NQH++B9E/wPt3wDfRkUvX4i7DGlZhfDIY+Tkc/9SrVYxpEWgyT7L398ftVpNVFQUTk5OefZFREQ88L1PPvkka9asKfAYtTrvDS9fX1+uXr1K8+bNAbhy5Uqe/RqNhieffJIVK1Zw/fp1nn/+eRmaYSEs8talp6cngYGB+jEdNjY2JCcn67+l3Wv58uV4eXnx7rvv4urqiq2tLc2aNSvhqIsmd0Xy/CjAtbhU03zQPxsgMQp8GwMqcCivG5P3dQfd9GJX/9IdJ1OMiYfQN8iXxx/1xkatIvfpsgrdrECPP+pN3yDTPa/y8fGhX79+TJw4kZs3bwIQHR3Nxo0bH/jeIUOGsGvXLtavX09WVhZZWVmEh4fz999/F/ieQYMG8cknnxAdHU1CQoL+OeDdcm9frl27lpEjRxb9lxMmZZGJLpdGo8He3p5+/foxdOjQ+25R5Prjjz/w8/OjR48eeHh40KRJE7Zt25bnmLNnz1KxYkWqVq3K+PHj9bcwzC13RfL8qIDK7k4F7DXC3c/kXtitm14sMwnqD4DyAbrpxb7tCiv66oYsbH1ZphgTRaJWq1gwqDFz+gfRONAdHzcHGge6M6d/EAsGNb6/c9VDWr58ORqNhmbNmuHm5kbbtm31z+EKU7lyZX799VeWLFlCpUqV8Pb2ZsKECSQmJhb4nnfffZegoCDq1q1Lw4YN6dmzJwAODg76Yzp06ICNjQ3VqlUjKCjo4X9BYRIlPjNK7969CQkJKXD/pUuXqFKliv51Wloa69evJz09ndGjR+f7ni5durB7927Wr19Pr169CAkJ4bnnnuP48ePUqFGDixcvotVqqVatGpcvX2bMmDE4OzuzefPmAuOYMWMG77//fp5txXGqNh29xmvr8r/VY6NWMad/0MPduiyo40nu9l5zQG0H++ZA3CXdPo/quj/nvkc6r5Q5MjPKg4WGhtKhQwfS09PzrCrQqVMnnnrqKSZOnGiSz5GZUR5eiSe6xMREMjMzC9zv4eFx371xgEcffZQlS5bQpk2b+/b169ePW7dusW/fPv22tm3bMmjQIMaPH3/f8adOnaJ+/fokJSXdd2+/IMV1AeQOsN1+MgbtXWOPABr4lWfT+NYP9y14z2xdKyy/BHVkha4XZ4cpkJMN/6yHfZ/BzX/vHONZE25fkCRXxkiiu19sbCwnT56kXbt2xMTE8Nxzz+Hr65unp2ZoaCjdu3fn6tWrlC9f3iSfK4nu4ZV4ZxQ3N7cHH5SPrKwszp07l2+iCwoKYteuXQaXlZtILeEvNPdWT/CxKP2K5Bone85EJ3HxRjI3kjPwditX9A8wtJemjS0EDYT6z8DpYN3sKrEn4dY5sHcGTYCub7gswCnKqJycHF555RXOnz+Pk5MTXbt25YsvvtDvf/zxxzl48CCff/65yZKcMA2LnNR569atBAQEULduXTIzM5k/fz4zZ87kxIkTVK9+/4rbFy5coEGDBqxdu5aePXuybds2Bg4cyPHjx6levTrbtm2jUaNGVKpUicjISEaNGoWtrW2ht1DvVdLfdGZuPcXS/ZfoUqci3wxtWvIrPB9eDiEvg205yErTbQtoBR2nQtW2hb1TWAFp0VkOadE9PIvsjHLz5k369++PRqMhICCAHTt2EBISok9y+/btw8XFRX989erV+fnnn3nzzTdxc3Pj3XffZf369frjd+/eTaNGjXBycqJly5ZUq1aNlStXmuV3M9Rr3WoR6OnEztOxJp826YGOrICQV6D35/B2FDT5b9DrlT9hRW9Y0QeuHCzZmIQQoogsskVniczxTSf0wi0GfXMQdyc7drzaHi8Xhwe/6WEV1PHk8DLY+go4e0HKDd226p10i8JWbvLgZ4CiVJEWneWQFt3Ds8gWndBpWd2TwS0CiEvNYkbwyZL50IKmGGs6QjcsoclIeHIJuFeFC7sgbIVuEdhd94wpkuEIQggLIS06A5nrm05yRjbd5+3lWnwaiwc34fF6Pg9+U0nIyYJja+CPTyDhvxkifBrAk4sh8rAMRyjlpEVnOaRF9/CkRWfhXBxs+egp3SKt0zb/Q3xqwUMzSpSNHTQeApOOQO95UE4D0cfhq1aw5SXoPF2SnBDCIkiiKwXaP1KB/k38uJGUwYdbjVyktbjZ2kPTkfD6v7pZVgBQYPdHsH0qpBixzp4QQhQDSXSlxLu96lLR1YH1YZHsORtr7nDud2wNJEb+N0H0f0MhDi6ELxrC3s90U4sJ6ydzpQoLJImulCjvZMfMfvUAmLrhBEkPWACzROWZS3OPrtOKNhuqtIWsVF1HlS8a63pu7vpIKkJr5lrJIudKzc7OtohnT3evcC5KjiS6UqTboz70CfIlKiGd2b+cMXc4OgWteN57Plw+AB3ehkefguTo/yrAZbD1JYurCIWJ5P7d353sinGu1CpVqjBr1iwaN26Mm5sb3bt3JypKN+5UpVLx5ZdfUq9ePZydnUlOTubChQv06dOHChUqEBgYyMyZM9FqtYBunt0uXbpQvnx5PDw8aN26NamputVD5s6dS0BAAK6urlSpUoWlS5cCujlx+/XrlycmjUbDnj179Pt79+7NuHHj8PDwYMqUKSiKwhdffEHt2rXRaDR06NCB06ct7JGElbHI9ehEwWb0qcv+czf48dAV/oq4TWJaFn4aR4a0rELfIF+Tzw7/QIWteJ67v/8yaDUJdr4Hl/bqtm+ZDLfOQ7cPZdLo0mJePUgveHb/PGwcdH/HW18BJQdsHeG3abqfBynnBq/8Y3BYS5cu5ZdffiEgIIBx48YxePBg/ZSAq1at4rfffsPT05OcnBw6d+7Myy+/zPr164mOjqZnz55UqlSJUaNG8c4771CjRg1++eUXAP7++29sbW35999/effddwkLC6N27drExMQQExNjcHzbt29n6dKlLFiwgMzMTL766iu+/fZbtmzZQtWqVVm0aBF9+vTh1KlT0sO1mEiLrpRxd7KniqczAOdikolJzCDsSjyvrTvGpNVhaPNZBaFYdZhScHJqMuzOYPHKjWFoMAzeAD66XqT8+QXMDpAkZ41s7EFlo0tyKhvd62Iybtw4ateujZOTE5988gm7d+8mMlK3wPCbb76Jr68vDg4OhISE4O7uzssvv4y9vT0BAQG89NJL+kmZ7ezsuH79OhEREdjZ2dGqVSvs7e2xsbFBURROnjxJWloa3t7eNGjQwOD46tWrx/Dhw7G1tcXJyYmFCxfywQcfULNmTWxtbZk8eTJpaWkcOnSoWM6PkBZdqRN8LIrj1xLybFOAHK3C9pMxBB+LerhlfYqTSgU1OkO1jrqVEjaPh/QEQKVbMSE9AcrJZLgWy4hWlr6VXrkpRIXpWu7F9EUmMPDOquXe3t44ODhw7do1AAICAvT7IiIi+Oeff9BoNPptWq0Wf39/AD799FNmzJhBly5dUKlUDB8+nOnTp1O9enVWrFjBl19+yYgRI2jRogWffPIJDRs2NCi+u2PIjWPw4MHY2Njot2VmZuqTszA9adGVMitDIwpstWm1CisPXi7hiIpArdZ1UtFmg1tlQIHQL+90WNHmmDtC8TDuvhU95vf7n9mZ2OXLd6752NhYMjIyqFxZ92Xv7iW//P39adKkCfHx8fqfxMRETp7UzTpUsWJFFi1axOXLl9myZQuLFy/Wr1Y+YMAAdu/eTUxMDEFBQQwZMgQAFxcX/XM8gJSUlPsWb7132TF/f3/WrVuXJ47U1FQGDRpkwrMi7iaJrpSJjE+joJuTCnAtLrWAvRbk7orw1VPQ9UNABak3dduXtIeI/eaNURRNYZ2TiinZLVmyhLNnz5KWlsZbb71Fu3bt8PPzu++43r17ExMTw6JFi0hPTycnJ4ezZ8/qO4789NNPXLlyBUVR0Gg02NjYYGtry9mzZ9mxYwdpaWnY29vj4uKCra3uZljjxo0JDQ3lzJkzpKenM3Xq1AeuNDJhwgSmT5/O2bNnAd0anZs3byYpKcm0J6YEBAcH07BhQ5ydnfH19WXx4sXmDilfkuhKGT+NIwX9M1IBld0NW0jWbPKrCFtP1g1JUKl1c2jGnIDlvWDtEPhligxHKE0K65zUe75uv4mNHDmSQYMG4e3tzbVr1/jxxx/zPc7FxYWdO3fy+++/U6VKFTw9PXnuueeIjo4G4MiRI7Rq1QoXFxdatmzJqFGj6Nu3L5mZmUybNg1vb288PT3ZtWsXy5cvB3Srib/44ou0atWKGjVqUL9+fVxdXQuNd+LEiQwfPpynnnoKNzc36tSpk2fx1tJi+/btjB8/nvnz5+tbxh06dDB3WPmSuS4NZClzwG06eo3X1h0jJ5/bl2oVzB3Q0HKf0cGDVzxPjAKvmrBjOiReA5WtrkNDj/9B8xfzHiudWIpNaZnrskqVKsyfP/++Lv7WxFLnumzWrBljxozhhRdeKNaYTEFadKVM3yBfHn/UGxu16r6Wnb2Nmg61KpglLoM9qJdmx7d1q5xPPKwbg2djByjwy5uwaTxotZLkhDCzlJQUjhw5wrVr13jkkUfw8fGhf//+XL9u+ha7KUiiK2XUahULBjVmTv8gGge64+PmQJNAdxoFaEjP1vLRNisZeGrvpEuKkw5D/f66beE/wmw/SXJCFCOVSqX/mTFjRr7HxMXFoSgKmzZtYseOHZw/fx4HBwcGDx5cssEaSIYXlEJqtYp+jSrnuUUZl5JJt/l7+elwJN0f9aFzHW8zRmhC5f3g6aXw2Avw/RN35sy8+S9kJIODS+HvF1YtIiLC3CFYHUNuXbq46P7dTZ48WT+84/3336dmzZqkpKTg7OxcrDEaS1p0VsLd2Z7Z/y3nM2XDCeJSLGQ5H1OJPQ3Z6XdWSAj9EhY+Bqe3gAU8OxWiLNFoNPeND8xlCX0Z7iWJzop0ruPNgKa65Xyml9SK5CXh7mdyr5z4bzgCus4qawfD6mchrhSMHyxFHtRFXpQ8S/s7eeGFF1iwYAHXrl0jLS2NDz74gM6dO+tbe5ZEEp2Vebd3XXzLl2PLsShCjlvmg2GjFDgc4QvdcAQnL/h3OyxsDvvmwq5ZMhzBBHJn7ZDZ9s0vNTUVlUqlH7tnKaZMmULnzp0JCgrC39+f1NRUVq5cae6w8iXDCwxkKcMLDLH/3E0Gf3sIdyc7fnulPRVcHcwdUtE9aDhC/BVQ28D+eZCTCS4+kBJ7f2cV6alptMjISDIzM6lcubLFtSbKitTUVKKjo/H09KRCheLvUV2a6jljSKIzUGm7AKZv/ofvQy/TpY433wxtYv0V1c3zEPIqXPrjvw0q6D4LWk6QJFdE2dnZXLp0iezsbHOHUmapVCo8PT3x8vIqkX/Dpa2eM5QkOgOVtgsgNTObHp/v4/KtVOb0D+LpJvdPiWR1FAVO/Ay/vg0pN3Tb3KvoWn2S5IpEURSLWbS0rMm9XVmSX1JLWz1nKEl0BiqNF8DhiNv0XxKKi4Mtv77cDl+No7lDKhlp8bpVzf/WLY5JufIwLhTKW/CMMUJYgNJYzxlCEp2BSusF8NG203y99yK1vF1wdrDlWnyaeRdqLSm5tytt7CA7A2zLQY9PoPFQ3XJBQoj7lNZ67kEk0RmotF4AqRnZNPtoJykZd5a+UaEbdP74o94sGNTY+pLd3c/kGgyANc/BBd2K01TrCH2/AE3+Y4CEKMtKaz33IDK8wMr9diqGtMy867vdu1CrVbm344mdIwzZCG1e1e2/uBsWtdTd1tRqzRmpEKKESKKzcitDIwqcOKTULNRqjIKWienyHvSaCwEtISsNQl6DL4Jg75z8y5Exd0JYDUl0Vs4qFmo1RmGrIzQbBSO3w5hd4F1P1xtz1wewbmTe1l1uq9C1UomELIQoXpLorFypX6i1OPg2hDG7dcsAqdRwcj0saKIbiydj7oSwOkZ1RklJSeGff/4hLi4Od3d36tWrZ3GzVBeX0vqQtrCFWm3UKub0D7LshVqLW/Q/sHogJERC7leCPp9LkhNlUmmt5x7EoBbdli1b6NatG+7u7nTt2pXRo0fTtWtXPDw86Nq1K5s3by7uOEURFbZQq5/Gkb5BvmaJy2L41IPJx8CtMrqbuQqc2gxJ0eaOTAhhIg9s0bVt2xZ7e3uGDh1K586d8fO7M8NGZGQku3bt4vvvvyczM5O9e/cWe8DmUpq/6Wi1CsHHolh58DLX4lKp6FaO87HJpGbmsGxEMzrWqmjuEM0r93ZlhdoQe0q3zdFDNwyhTh+zhiZESSrN9VxhHpjo/vzzT1q1avXAgkJDQ2nZsqXJArM01nYBbP8nmrE/HKGCqwPbX2qLp0spnvj5Ydz7TO6vpbDtdcjtwtNoMDw+GxxczRikECXD2uq5XDJg3EDWeAG89fNx1h6+Sre63iwZUgYmfr5XQR1PcrfbO0NGkm6+zKe+Af/HzBOnECXEGus5KEKiCw8PJywsjOTk5DzbJ0+ebNLALI01XgApGdn0/EI38fPsp+rz7GNlbLaQBy0BdOsC3L4AZ7aCykY3Bq/+M9B0RP7HJ13XDW8QopSyxnoOjEx07777LnPmzKFBgwY4Od3plq5Sqdi1a1exBGgprPUCCLsSR//FodjbqNn2UluqepWNXrQGUxQ4+gP88hZkpei2dZoG7V6/c4wMSRBWwlrrOaMSnaenJ/v376dOnTrFGZNFstYLAGD+zn+Zv/McQf4afh7bEjsbGV55n9sXYcMLEPm37nWDZ+HJxRD2vSQ5YTWstZ4zKtHVrFmTEydOUK5cueKMySJZ6wUAkJ2j5ZnFoYRfjeelzjV5pesj5g7JMuVkw745sOdjQIFyGshIlCQnrIa11nNGfXX/7LPPmDBhAufPnycxMTHPjyi9bG3UzB/YECd7G77cfZ4jl+PMHZJlsrGFDm/BqB1g4wDp8aC21Q1LEEJYLKNadH/88QeDBw8mKurOjPeKoqBSqcjJySnknaWftX7Tudvav6/w1voTBHo6ETK5LS4OtuYOyTLlPpMrp4G024AKusyAVpNBLbd9RellCfVceno6f/zxB+Hh4fpZuBo2bEi7du1wdCza4tFGJbrq1aszePBgnn322TydUQACAwOLFEBpYQkXQHFTFIUXVx7ht1MxNK/qQXaOlsiyslCroe7ueNJ4KGwaB8dW6/bV6ApPLgFnT3NGKESRmbOeu337NjNnzmTZsmV4enpSt25d3NzcSExM5NSpU9y6dYvhw4fz7rvv4ulp3L8xoxKdRqMhLi6u7I23omwkOoCbyRm0nr2LjOw7s/lb/UKthiqod+Xu2fDHx7o/u/rCM99BoPVOniCslznruRo1ajBs2DCGDh2ab8PpypUrrFixgpUrV/Lvv/8aVbZR91meffZZgoODjfqAoggLC6NJkyZ4eHig0Who1apVodOLjR07FhcXF/2Pk5MTKpWKsLAw/TGbNm2iZs2aODk50aZNG86cOVPsv0dptP/cTbJy8i5IatULtRqjoLXuOk6BHp/qlv5JioLlvXSdVmRhVyEMFh4ezrRp0wq8OxgQEMC0adPy1OuGMqpF16dPH3bs2EGTJk3w9vbOs2/Dhg1Gf3hBbt26RXJyMgEBugHMGzduZNiwYcTGxhp0j3bOnDl8/fXXnD17FoCzZ8/SuHFj1q5dS5cuXfjoo49Yu3YtJ0+exNbWsOdQZaVF9/SiA4Rdic93DTsV0DjQnfXjHjwlXJl19Efdoq7ZaeBeFZqNhlYT7z9OBpgLC2St9ZxRvQ2aNm1K06ZNiysWPU9PT/09WK1Wi42NDcnJyURHR1O1atUHvv/bb79l5MiR+tc//PADHTt2pHfv3gBMmzaNBQsWsG/fPjp27Fg8v0QpVeYWajW1Rs9D5cawbjjcOAO/vQPJMdDtwzvH3H0LVAiRr9DQUPbs2cPNmzfzJN+5c+caXZZRie69994z+gMehkajITk5mZycHIYOHWpQkgsNDeXcuXMMHz5cv+348eM0bNhQ/9rOzo66dety/PhxSXT38NM4EpuYUWCLrkwu1GqsinV0q5hvexPCf4A/v9BNJzZwpW6WFRlgLkShFi5cyBtvvEH37t355Zdf6NGjB7/99htPPPFEkcozui90aGgoY8eOpXfv3owdO5bQ0FCj3t+7d29UKlWBPxEREfpj4+PjSUpKYuXKlbRt29ag8pcuXUrv3r3z3FpNTk5Go9HkOU6j0ZCUlFRgOTNmzMgTV1kxpGWVAjubqFUqhrSw7t61JmPvDP0W6nph2tjD2RD4XxVJckIYYP78+fzyyy9s3LgRR0dHNm7cyLp163BwKNoqK0YlujVr1tCtWzcURaFt27aoVCoef/xxVq9ebXAZq1at4saNGwX+5D6Xy+Xo6MjgwYOZN28e+/fvL7Ts5ORkfvrpJ0aNGpVnu4uLCwkJCXm2JSQk4Opa8NIrM2bMQFEU/U9ZUdhCrRXdHOjToJJZ4iq1gp6FsQfA1lE3i4raFnzqmzsqISxaTEwM7du3B0CtVqMoCj169ChyZ0ijbl3OnDmTkJAQ2rVrp9/23HPPMXbsWAYNGmRQGW5ubsZF+J+srCzOnTtHmzZtCjxmzZo1uLm50aNHjzzbGzRoQHh4eJ6yTp06Rf36UuHcS61WsWBQ4zwLtfqUd+TyrRSuJ6SzPuwaA5r5mzvM0uVKKORk6BZzTbsNS7tCn/nQeIi5IxPCIvn4+BAVFYWvry9VqlRhz549VKhQAXVRJ2RQjKDRaJTs7Ow827KzsxWNRmNMMQ+0ZcsW5dixY0pWVpaSkpKizJo1S3F0dFTOnz9f6PtatGihTJ069b7tZ86cUZycnJSQkBAlPT1dee+995SaNWsqWVlZBsdk5KmyOocjbinV3g5Rar27TTkXk2jucEqPw8sVZYZG93+tVlF+GqEo77npfoInK0pWurkjFELPUuq5uXPnKhs2bFAURVF++OEHxdbWVrG1tVXee++9IpVn1G/VunVrZeHChXm2ffXVV0qrVq2K9OEFWbZsmfLII48ozs7Oiqenp9KhQwdl165d+v179+5VnJ2d87zn5MmTikqlUi5cuJBvmRs2bFBq1KihlCtXTmnVqpVy+vRpo2KylAvAnBb8/q8S+NZWpfu8P5S0zOwHv6GsuzvJ3e23aXeS3dcdFSU+0jzxCXEPS6nnbt++nef1lStXlFOnThW5PKPG0R0+fJgePXpQsWJFqlSpQkREBLGxsfzyyy8lMuzAnKx1fIkxcrQKQ749xJ8XbjGkRSAf9qtn7pAsW2ELux74Ag4thsRr4OQF/ZdDVcM6XAlRXMxdz4WFhfHEE08QFRVFYGAgwcHB1Kv38PWM0SuMJyQkEBISQmRkJH5+fvTs2fO+Ho3WyNwXgKWISUynx+f7uJ2SyeLBTXi8no+5Qyq9sjPh17fh76W6Fcy7fgAtJ0AZ6uUrLIu567lOnTrRsGFDRo0axTfffMO///7Ltm3bHrpcoxLd7NmzmTLl/pkcPvnkE958882HDsaSmfsCsCS7z8QyYvnfuJWz5ZeX21FZU7QZxcV/wlfB1lcgOx18G8GwreDgkvcYmUlFlABz13NeXl5ERUVhb29PamoqNWrUyLNaTlEZlehyZ5K+l6enJ7du3XroYCyZuS8ASzNz6ymW7r9E00B31rzQAltZlfzhXD8G3/fT9cp08YER28Czum5fQZNJC2Fi5q7n7s0xHh4e3L59+6HLNWh4wfHjxwHddFwnTpzIcyIuXLhQ5DWCROn15uO1OXTpNocvx/HF7+d4tVstc4dUulUKgklHYFkP3dRhX7WCZ5ZByg1JcsJipaWlUb9+fW7evEl8fPxDl5eZmckXX3yhf52enp7nNcDkyZONLtegFp1ardbPDnL34SqVCh8fHz788MM8c0taI3N/07FEETdT6PXFPlIycxjfsTqHLtyS9eseljYHVg+Eczv+26CCPp9LkhMlwth67o033iAsLIwjR46YJNF16NCh0JmoVCoVu3btMrpco25dNm/enEOHDhn9IdZAEl3+NoRF8upPxwDdXJgKsn6dSXzRGG5f0P253jPwxJdgJ3dORPEypp47cuQIw4cPZ86cOQwYMMAkia64GPVgpawmOVEwterOVGHKXf+X9esewpEVEHcJKtbVvf7nZ90tzUQ5l6L43T3H74wZM/I9Jjs7mzFjxrBw4ULs7e1LNsAiMGoKMIDVq1fnu3SCKdejE6XHytCIAvdptQorD16mX6PKJRdQaXdvx5M/v9Qt9RN1FL7uCM+uAr8m5o5SWDFDWnSffvopjRo1ol27duzZs8dkn/3nn3+yZcsWPv74YwB8fX1JT0/X79+6dSutWhm/HqZRLbrp06fz6quv4uTkxC+//EJgYCAHDhzA31/mPiyrZP06E8qvd2WridBrLqCC5Ghdy+74T2YMUpR158+fZ/HixXz66acmL3v+/Pl5Jh9JSUlh48aNbNy4kbfeeot58+YVqVyjntFVqVKF4OBgGjRogEajIT4+nkOHDvHxxx+zadOmIgVQWsgzuvzJiuQmVNhMKkdWwOlguLALFC20eQU6TYeiTnIrRD4MqeeWL1/O2LFjcXHRjfXMysoiKSkJDw8PQkJCaN68eZE/v2rVqvzzzz84OzsD4O7uTlxcHKDr4fnoo49y8eJFo8st8ji6ChUqcP36dWxtbfVJz5pJosvfpqPXeG3dMXK0958bG7WKOf2D5NalKZ3/HX4eAekJ4FkDHnsBmr94/3EywFwUgSH1XGpqap6xbaGhoYwePZqTJ09SsWLFh3pmd+84ul27dtGpUyf9a1dX10LXES2IUV8HAwMDOX/+PACPPPIIa9asYdu2bfrsK8qewtavq17Bmb5BvmaJy2rV6Ayjd+mS3K3z8MubsHdO3mNyb4G6ytqBwvScnJzw8/PT/1SoUAGVSoWfn99Dd0xxc3Pj8uXL+td3J7mIiIhC1xAtjFEtuh9//BEvLy+6d+/Or7/+ytNPP01GRgYLFy7khRdeKFIApYW06Aqm1Sp51q/zdHHgXEwSmTkKP4xqTpuaXuYO0fqkxcPPI+HC77rXLSdB95kyi4p4KOau50aNGoW9vT1fffXVffsmTJhASkoKy5cvN7pcgxKdoij5DuLLysoiIyNDf6/Wmpn7AihtNodf46U14Xg62xMyuS0+5cuZOyTrk5MNO6bDwYW615pASLgqSU4UmbnruatXr9K0aVMaNWrEs88+i6+vL1FRUaxbt47Dhw/z999/ExAQYHS5Bt26rFixIsOGDWP9+vWkpKTot9vZ2ZWJJCeM90TDygxuEcCtlEwmrQ4jK0dr7pCsj40tPP4R9P0SUEH8ZXCuAI0GmzsyIYrE39+fv//+Gx8fH6ZOnUrv3r2ZOnUqXl5eHDp0qEhJDgxs0UVFRbFlyxa2bNnCgQMHaN68OX379qVv3774+fkV6YNLG3N/0ymNMrJz6L84lOORCbzQrhpTe9Yxd0jWKfd2JSpQcsCnAYzcDvby7FwYx1rrOaPXo0tNTWXHjh1s2bKFkJAQfHx89EmvSRPrHchqrRdAcbt6O5XeC/aTkJbFkiFN6P6orF9nUnc/k6vSBr7tBqk3obw/jN4JrnK+heHMWc/dunULT09Pkx13N6MT3d0UReGvv/5iy5YtBAcH8+yzzzJ16tSiFmfRJNEV3e+nYxi14jCu5WzZOqkNgZ7S0jCJ/DqepN6GpV1082SWc4cRIeD9qDmjFKWIOeu5GjVq8OSTTzJixAjq1q173/7Tp0/z3XffsWnTJs6dO2dU2QYluuHDhzNmzBhat25d6HFZWVnY2dkZFUBpIYnu4fxv+xm+2nOBR33dWD+uFeXsbMwdUulX0ADz7AxY1hOuHQZ7VxiwHGp0MUuIonQxZz2XlJTEnDlz+Prrr9FqtdSuXVs/ru7s2bMAvPjii7z22mtGDzMwKNE9//zzbNq0CX9/f0aNGsXQoUPx9vYu2m9TSkmiezjZOVqeX3qIQ5duM+ixAD5+qr65Q7JuigJ//A/2fAwqG+j1GTS17qW0xMOzhHouJyeHv/76i/DwcOLi4nB3d6dhw4Y89thj2NgU7QuywbcuExMTWbVqFcuWLSM8PJwePXowevRoevbsiboMTENkCRdAaRebmE7PL/ZzMzmD55sHcOZ6oqxfV9yOrYXgiZCTCS0nQtcPZdowUSBrreeK9Izu5MmTfPfdd/z444+o1WqGDx/ORx99VBzxWQxrvQBK2oHzN3l+ad7lnmT9umIWcQDWPg9pcf/1yPwV7J3yHiNThgmst54r0le7Rx99lDlz5rBq1SrKlSvH//73P1PHJazUjaQM7p17QNavK2ZVWsPo38HJC6KPw8LmkBRzZ79MGSasnNGJ7vr163z88cfUqlWLJ554wuTrEQnrtjI0goLW9cldv04UA8/qMPFv8KgOCVdgUXOIOSVThokywaCFV7Oysti8eTPLli1jx44dNGrUiNdee41BgwYVeZJNUTbJ+nVm5OQB40Pv9Mhc3AZQJMkJq2dQoqtUqRJqtZrnn3+e//3vf9SrV6+44xJWyk/jSGxiRoHr11V2d8pnjzAZWwfdQPK5dXTP5FDJDCrCIu3evZtVq1YRHR3Nli1bOHz4MElJSXTs2NHosgy6dfnVV19x7do15s2bJ0lOPJQhLasU2NlEpYIhLQJLOKIyKOx7SI7RTQKNAutHwYEvdEMShLAAS5cuZciQIXh7e7N3715AN7fy9OnTi1Se0b0uTZllSxNr7Y1U0rRahUmrw9h+MgatVsnTsnNxsCF0SmdcHa1z0gGLcO8zuW1vwF9f6/Y1HwvdPwK1DOYvqyylnnvkkUfYtGkTdevW1a8ynpWVRaVKlbh586bR5RnVGcXUWVaUPWq1igWDGjOnfxCNA93xcXOgSYCGWt4uJGfk8NaG4xbxD80q5dfxpOen0OY13Z8PLYZ1wyErzUwBCqFz69Yt/TRguUvEqVSqfJeLM4RRLTpTZ9nSxFK+6Vir+NRM+n55gCu3U3mjey0mdKxh7pCsT0FThgHs/Qz+XADp8eDfAgat1nVeEWWKpdRznTp1YsKECTz99NN4eHhw+/ZtNmzYwJIlS/j111+NLs+oROfp6cmtW7cA9B+enZ1NpUqVuHHjhtEfXppYygVgzc5EJ/LUoj9Jy8rhu+HN6FirorlDKluSb8CqARAVBl6PwPM/g7s8My1LLKWeCwsLo1u3brRv356QkBCeeuopdu3axa+//kpQUJDR5Rl16zIoKIj169fn2RYcHEzjxo2N/mAh7lXbx41PnwlCUWDy6qNcupny4DcJ03GpAMO3wiOPw81/4auWsGtW/sceWaFrIQpRDBo3bszJkydp2bIlo0ePJigoiLCwsCIlOTCyRWfqLFuaWMo3nbIgd6WDmhVd2DihNS4OBo2CEaaSkw0hr0LYCt3rFuPh8Y/v7JdB5lbLEuq57Oxs/P39uXTpEuXKlTNJmUa16EydZYXIz+vdatH+kQqci03m9Z+Omf0fXpljYwt9PoeO7+heH1wEmybo/ixJThQzW1tbXFxcyMzMNFmZD7XwalliCd90ypKE1Cz6LtzP5VupvN7tESZ2qmnukMqmoz/A5omAAm6VdYPMJclZLUup55YvX86WLVuYMWMG/v7+eVbIcXNzM7q8Bya6xYsXM3r0aGxtC759lJOTw9KlS3nxxReNDqC0sJQLoCw5G53Ek4sOkJqZw5i2VQm7HCfL+pjD+Z3wY39QtODqA6+e4b6ZuYVVsJR67u7EljukQFEUVCoVOTk5Rpf3wIcfx44do3r16gwcOJBOnTrlWfX1zJkz7N69mzVr1tCzZ0+jP1yIwtTyceXTZxowYdVRvtl3Sb89NjGD8Mhj7DgVLcv6lISEa7r/q2wgKRpW9IahW2RdO1FsLl269OCDjGDQrcvz58+zZMkSgoODOXfunD7D1qxZk969e/Piiy9Ss6Z131qylG86Zc2mo9d4ZW14vnNj2qhVzOkfRL9GlUs8rjLj7mdylRvDd49DZjJUbgYjt+ue5wmrYa31nNHP6NLS0vTLmzs6OhZXXBbHWi8AS/f0ogOEXYkvcBLoxoHurB/XqqTDKhvy63hy8zws7QLp/y3iOnqnbqJoYRUspZ4bOXJkgfu+++47o8sz+t6Do6Mjvr6+ZSrJCfORZX3MKL+OJ141YOy+O4u4rhoImTLeUZhW+fLl8/ykpqby888/Y2dXtHlw5b6DsGiyrI8ZdZiS/3aNP4z7E1Y+CRd3w8qn4Lm14Kgp0fCE9Zo3b95923bv3s0333xTpPLkabKwaLKsj4Vy9dbNolK5KVw9CCv6QIp1z3crzKtDhw6EhIQU6b2S6IRF6xvky+OPemOjVnFvulOrVNT1NX5MjTARJw8YugmqtNXdxlz4GOz/PP9jZcowYYTExMQ8P9HR0XzyySf4+PgUqTyjOqP8/fffNGvWrEgfVNpZykPaskirVQg+FsXKg5e5FpdKZXcn3J3s2Hk6lkBPJzaOb42Hs725wyy7stJ0y/v8u133utM0aPf6nf0ym0qpYSn1nFqtzrMkj6IoBAYG8t133xVp7VOjEp2HhweVK1dmxIgRDBkyhAoVKhj9gYYICwtjzJgxXLp0Ca1WS926dZk9ezbt2rXL9/ixY8fyww8/6F9rtVrS0tI4cuQIjRs3Zs+ePXTs2BFnZ2f9McOHD+fLL780OCZLuQCEjlarMHF1GNtORPNYFQ9Wjn4MB1tZMNRscrJg44vwz3+Tvnd4W/eMT5JcqWIp9dzly5fzvHZxccHT07PI5RmV6DIyMtiwYQMrVqxg7969dO/enZEjR9KrV688I9kf1q1bt0hOTiYgIACAjRs3MmzYMGJjYw3q7Tlnzhy+/vprzp49C8CePXvo168f8fHxRY7JUi4AcUdaZg4Dvw7leGQCTzf247P+DYq8MKMwAW2OLqmFfa97XbEO3DgrSa4UsdR67ujRo9ja2lK/fv0ivd+o7OTg4MCgQYPYvn0758+fp3nz5rz22mtUrlyZN998kwsXLhQpiHt5enoSGBioP+k2NjYkJycTHR1t0Pu//fbbQsdhCOvgaG/D0qFNqVS+HOvDIvnqD9Ncf6KI1DbQ5wtoOVH3OvY0VKgtSU4Y7YknnmD//v0ALFy4kFatWtGyZUsWL15cpPKK3Aw7f/48Z8+eJSYmhrp163Lr1i2aNGnC7Nmme+Cs0Wiwt7enX79+DB06lKpVqz7wPaGhoZw7d47hw4fn2Z6cnIyvry9+fn48//zzXLt2rdByZsyYoV+6XVoJlquiWzmWDmuKk70Nn2w/yy8nrps7pLJNpdIt2prbdSj2FOz5n1lDEqaXkZHBmDFjqFq1Kq6urtSuXbtIA7kLEhoaSvPmzQFdotu5cycHDx5kzpw5RStQMUJERITy/vvvK9WqVVP8/f2VadOmKRcvXtTvP3/+vOLm5lZoGb169VLQjfXN9+fSpUt5jk9NTVVWrlypfPPNNwbFOHLkSKVfv355tl2/fl05ceKEkp2drVy/fl0ZNGiQ0qhRIyUnJ8ewX1xRFCNPlShhO05GK1WmbFVqvbtNOXY1ztzhlF2HlyvKDI3u/9veVJT33HQ/ez4xd2TCAIbWc8nJycq0adOU8+fPK1qtVgkNDVU0Go3y66+/miSO3DwSGRmpeHt767e7uroWqTyjntE5OjrSp08fRo4cSffu3fNt6YwdO7bQ5mViYmKh6wx5eHjk+7zv0UcfZcmSJbRp06bA9yYnJ1OpUiVWr15N7969Cz2ufPnynDx5ktq1axd43N0s9d61uGPpvovMDDlNBRd7JnSswZZjUbLaQUm6t+OJosDWV+DIMt3+TtOh3WvmjFA8wMPUc0899RT16tXjgw8+eOg4WrZsSd++fbl8+TJpaWmsWLGC2NhYgoKCuH7d+Ls2Rs2McvXqVby8vAo95kH3UIuylhBAVlYW586dKzTRrVmzBjc3N3r06FFoWXIr0jqNalOV87FJrPk7khlbTqFCd5tAVjsoIfdOGaZSQa+5kJMJ4T/Cgc+hwQDdzCrCYt1dP7733nvMmDHjge9JT0/nr7/+4rnnnjNJDIsWLWLixInY29vrb4n++uuvdOvWrUjlGdWi27t3b77bHRwcCAgIoFKlSkUK4l5bt24lICCAunXrkpmZyfz585k5cyYnTpygevXqBb6vZcuWdOrUiVmzZuXZvnv3bqpUqUKVKlW4ffs2r7zyCuHh4Rw9ehQbG8O6pEuLrnT4+chVXl93PN99stqBmWhzYMML8M/P4F4VRvwCbqapK4RpFaWeUxSFIUOGcO3aNX7//XeT9sA3FaMieuKJJ+jcuTMdOnTg8ccfp0OHDnTu3JkuXbrg5+dHu3btuHr16kMHdfPmTfr3749GoyEgIIAdO3YQEhKiT3L79u3DxcUlz3tOnTrFoUOHGDVq1H3lHT16lHbt2uHi4kK9evXIzs5m69atBic5UXqsPnTlvhlUcmm1CisPXi5gryg2aht4cgnU6Qtxl+D7vpAca+6ohAkoisL48eM5e/YsmzZtMlmS27ZtG+fPnwcgIiKCJ554gqeffprIyMgilWdUi27RokWEhYUxe/ZsvLy8uHHjBu+88w6NGjXiySef5KWXXiItLY3g4OAiBWPJpEVXOjT/aCcxiRkF7vdxc+Dg1C4lGJHQy86En4bCv79AxbowbCs4F30QsDA9Y+o5RVGYMGECBw8e5Pfff8fd3d1kcdSuXZudO3fi5+dH//79sbe3x9nZmevXr7NlyxajyzMq0fn5+XHhwgUcHO6sP5WWlkbNmjWJjIwkLi6OmjVrcvOm9U3uKomudJD16yxcdgasHgQXfgef+jBsCziaroIUD8eYem7ChAns37+fXbt2PdSsJfnRaDTEx8eTnZ1NhQoVuHLlCg4ODvj6+hYpvxjVzszOzr6v6RgVFUVWVhYArq6uZGdnGx2EEKYiqx1YOFsHePbH/yaCPgGL20F64v3HySTQFu3y5cssWrSIs2fPEhgYiIuLCy4uLowdO9Yk5Ts6OhITE8OePXuoXbs2rq6uqFQqfa4xllG9LkePHk23bt2YMGEC/v7+XL16lUWLFjFmzBgAQkJCqFOnTpECEcIU+gb5suNUNNtPxqDVKve17AI8Zf06s7Nz1K1ft7gN3L4IS9rrFnN1+O+5+93DFIRFCgwMLNY7XEOHDqVZs2ZkZGTw/vvvA3D48GGqVatWpPKMunWpKArfffcdq1evJioqCl9fXwYNGsTIkSNRqVTk5OSgKAq2tta3nqvcuiw98lvtwLd8ObYcv46Hsz3rx7WiqpfzgwsSxSs9UZfs4i+DZ014cS+cWCeTQJuRJdVzO3bswM7Ojg4dOgC6RJeYmEinTp2MLsvgRJednc2gQYNYuXIl5cqVM/qDSjtLugCE8RRF4YOtp1h2IAJ/D0fWj2tFRdeydx1bnLR4WNwaEiLBwQ0ykyXJmZGl1XPXrl3j6tWrtGjR4qHKMfgZna2tLfv377fK1pqwfiqVimm96tKrQSWu3k5jxLK/SUov2v1+YUKOGnhxH9g6QkYi2LtC0LPmjkqY2fXr12nfvj2BgYF06aLrJf3TTz/pH5MZy6jOKC+++CKffPJJkT5ICHNTq1XMHRBEy2qenIxKZOwPR8jM1po7LHF6C+RkgG05yEiAbzrphiKIMmvcuHE0b96c5ORk7OzsAOjcuTO///57kcoz6hldo0aN+Oeff9BoNFSuXDnP4MCwsLAiBVBaWFqTXhRdYnoWA5cc5PT1RPoG+TJ/YEOZFsxc7u548sjjumd2KbHg0wDG7AIbO3NHWKZYSj1XoUIFrl+/jq2tLR4eHty+fRuA8uXLk5CQYHR5Rt2HfPnll43+ACEsjVs5O1aMaMaTi/4k+FgUiWlZJKVnyQTQJS2/1cdf3AuL20L0cVjaBUb/DjbyuKSsKV++PLdv36ZixYr6bZGRkXh7exepPKNadGWZpXzTEaZzLjaJnp/vIyvnzt+rCt0tzscf9ZYJoIvbntngWun+jicJ12BJW0i9BfX766YPU8t0fSXBUuq56dOnExoayty5c2nXrh2HDh3i1VdfpVWrVkydOtXo8oyemGzZsmV06dKFBg0aAPDHH3/w008/Gf3BQpjbyWuJ5Gjz/qNWgBytwvaTMQQfizJPYGVFhyn5964sXxle2APlA3TDDTZPBK08Sy1Lpk+fTqNGjWjVqhUJCQk0adKEOnXq8OabbxapPKMS3axZs5g3bx7PPvssV65cAaBSpUp8+umnRfpwIcxpZWgEBX15lQmgzUwTAMOCwa0yHFsFWyZLsisjsrOzWbJkCR988AFJSUnExsaSmJjIp59+WuRe/0bduqxatSr79u3Dz88Pd3d34uLi0Gq1eHl56R8WWitLadIL05EJoEuBWxdgWU9IjoamI3Xr28l6ksXGUuq53LkuTcWoFl1KSop+zbncxfmysrLyTPIsRGnhp3EscEkfFVDZXaYLMzvP6jB8K9g5w+Hv4Je3uK8ZLvNiWp3u3buzc+dOk5VnVDuwRYsWLFq0iEmTJum3fffdd7Ru3dpkAQlRUoa0rEJ45LH7ntOB7lldu5peJR+UuJ9XTWj9Euz5CP5aohty0G2mrmUn82JaJVdXV/r160e3bt0ICAjIM5Rt7ty5Rpdn1K3Lixcv0rlzZzw8PDhx4gRNmzYlJiaGnTt3UrVqVaM/vDSxlCa9MB2tVmHS6rA8E0Cr0P1HUcDVwYY1L7bkUd/yZo5UALD7Y/jjv5Zb65fBo5rMi2lillLPjRgxosB9y5YtM7o8o4cXpKWlERISQkREBP7+/vTu3RtnZ+ufINdSLgBhWvlNAD24eQBnY5JY/MdFPJztWftCC2p6u5o7VAGwaxbszZ2dSQV9PpckZ0KWUs/FxcWZdCFXGUdnIEu5AETJUBSF97ecYvmfEVRwdeCnF1vKigeWYlELiD2t+3P7t6Cj8eOqRP7MXc+FhYXxxBNPEBUVRZUqVQgODubRRx996HKNSnRJSUnMmzePI0eOkJSUlGffrl27HjoYS2buC0CUPEVRmLrxH1b/dQXf8uVY+2JL/D2kg4pZ5T6Tq1DrTrLrMBU6vGXWsKyFueu5Tp060bBhQ0aNGsU333zDv//+y7Zt2x66XKMSXb9+/bh69SpPP/30fbcrX3rppYcOxpKZ+wIQ5qHVKrz+8zE2hF3D38ORdS+2wqe8LO9jFvdOGbbzfdj/X8eETu9CuzfMGZ1VMHc95+XlRVRUFPb29qSmplKjRg2ioh5+4gajEp1Go+HKlSu4ubk99AeXNua+AIT5ZOdoeWltOCHHr1PVy4kRrauy+eg1mRuzJOU3LybAjhlwYJ7uz53fg7avmiE462Hues7NzY3ExET967sndH4YRg0v8Pf3JytL1vASZYutjZr5AxuSnpnD72dimb75JCp0QxBiEzMIjzzGjlPRMjdmcUq6nn/vyq4zAC2ELoLf3weVGtq8XPLxCZPIzMzkiy++0L9OT0/P8xpg8uTJRpdrVItuyZIlrFu3jrfeeuu+WaRz5760Vub+piPMb93hq7zx8/F899moVczpH0S/RpVLOCoBwOU/4YdnICtFN8au1aQHv0fcx9z1XIcOHfSTkeRHpVIVqT+IUYnu7kF79354Tk6O0R9empj7AhDm9/SiA4RdiSe/q0AFNA50Z/24ViUdlsgVsR9+7A9ZqdD9I2g5wdwRlTrWWs8ZNQWYVqvN98fak5wQAJHxafkmOdDdxrwWl1qS4Yh7VWkDz60FW0f4dSocXGzuiISFMHqZnvzExcWZohghLJrMjVkKVG0Hz60BtS1sfwv++ub+Y2RuzDLHoER39yqvoJtw826BgYGmi0gICzWkZZUCO5soQJ8GviUbkMhftQ7QbIzuz9teh7+X3tmX23vTtZI5IhNmYlCvy7S0tDyv//777zyvrfGerhD36hvky45T0ffPjYku0X27/yKd61SUQeWWoMd/LbZDX0HIa6Cy0fXIlLkxyySDOqM8aGzDvfutkbU+pBXGyW9uzEGP+bPrdCzb/ommUvly/Di6OdUquJg7VAGw7S34K/dZncyN+SDWWs9JojOQtV4AwjSyc7S8tf4E68Mi8XJx4MfRzanlIxNBW4QFTeHWOd2f+3whia4Q1lrPGXTr8kGD+GQQuSjrbG3UfPpMAxzt1fxw8ArPfh3KylHNqVdZlvgxqyMr4PYF8KwBt87Dlv8GG0uyK1MMatE9aBAfwO7du00WlCWy1m86wrQUReGjbaf5Zt8lXMvZsnzEYzQJNN1yI8II904btu0N+Otr3T5p2eXLWus5WabHQNZ6AQjTUxSFeTvP8cXv53C0UzOiTVUOXbglc2OWpILmxpRkVyhrreck0RnIWi8AUXwW7T7PJ7+ezbNNBajVKh5/1FvmxixOe2brhhDkl8i2vaEbcqBooe8CaDy05OOzUNZaz5lkwLgQ4n6+GkfuveOvADlahe0nYwg+9vDLj4gCdJhScGut56fw7GpQ20HwJAj7vmRjEyVOEp0QxWRlaAQFzRmm1SqsPHi5ROMRd6n1OAz8QZJdGSGJTohiInNjWrhaj8PAlf8lu8kQttLcEYliIolOiGJS2NyYoLu1KcysVg9dslOpIHhi/slO5sbMV1ZWFhMnTsTd3R0PDw8mTZpEdna2ucPKlyQ6IYpJYXNjAsSnZRGfmlmCEYl81eoBzUbr/nxvspO5MQs0c+ZM9u/fz6lTpzh58iT79u3jo48+MndY+ZJelway1t5IovhotQqTVofdNzemWqXCxcGGhPRsqldwZvmIx2R+TEtw99CDvl/qemWWsbkxjann/P39mTdvHs888wwA69at4/XXX+fyZct79iyJzkCS6ERR5Dc35pAWgXSuXZGX14bz+5lYvFzs+W54Mxr4acwdrrg72ZXBuTENrefi4uLw8PDg3Llz1KhRA4Bz587xyCOPEB8fT/nyljUjkNy6FKIYqdUq+jWqzPpxrTg4tQvrx7WiX6PKuDrasWRIEwa3COBmciYDlxxk56kYc4cren6qmy4MAEXXqitjVCqV/mfGjBn5HpOcnAyARqPRb8v9c1JSUjFHaDyD5roUQpierY2aD5+oh5+7E7N/OcMLKw/zXp+6lHe0Z2VohMykYg5HVsDti3fmxtz6MqBA05HmjqzEGNKic3HRrc6RkJCAl5eX/s8Arq6WN5m5RbbowsLCaNKkCR4eHmg0Glq1asXevXsLfc+sWbMIDAzEzc2NRo0a8dtvv+XZv2nTJmrWrImTkxNt2rThzJkzxfkrCGEQlUrF2PbVWTCoETYqFe8Fn+KVteGEXYknJjGDsCvxvLbuGJNWh6HVyq3zYnX3tGGTjsBjY3Xbt76Sd/FWgbu7O35+foSHh+u3hYeH4+/vb3G3LcFCE11gYCAbNmzg1q1bxMXF8frrr9OrV6/7FoDNtWnTJj777DO2bt1KQkICr776Kk8++aR+KaGzZ8/y/PPPM2/ePG7fvk2nTp144oknLLYrrCh7+gT5MrZDdUA3xi43pclMKiUkv7kxe/4Pmo/T/TnkNfjrG3NFZ5FGjBjBrFmziI6OJjo6mo8++ojRo0ebO6x8WWSi8/T0JDAwUP9g1MbGhuTkZKKjo/M9/uLFizRr1oz69eujUqkYMmQIWVlZXLx4EYAffviBjh070rt3b8qVK8e0adOIjY1l3759JflrCVGoP8/fLHCfzKRSzJKu59+7ssdsXbJT28K21+Hg4nzfXhZNmzaNli1bUqdOHerUqUPr1q2ZOnWqucPKl0UmulwajQZ7e3v69evH0KFDqVq1ar7HDRw4kOjoaI4ePUpOTg7Lli3Dz8+PevXqAXD8+HEaNmyoP97Ozo66dety/Pjxkvg1hDBIZHz+dyxAZlIpdoXNjdljNjy/DmzLwfa3IHRRycZmoezs7Fi4cCFxcXHExcWxYMECbG0ts9tHiSe63r175+nVc+9PRESE/tj4+HiSkpJYuXIlbdu2LbDMihUr0qtXL5o2bYqDgwMvv/wy33zzDeXKlQN0PYTu7h0EuiRaWO+gGTNm5IlLiOL2oJlUKrvLWDuzqd4JnlsLto7w69sQutDcEQkjlHiiW7VqFTdu3CjwJyAgIM/xjo6ODB48mHnz5rF///58y/zggw/Ytm0b//77L5mZmWzevJmBAwfqH5S6uLjoewTlSkhIKLR30IwZM1AURf8jRHF70EwqKiAtM6fkAhJ5VetwV7KbCj8V0AKUKcMsToknOjc3N7y8vAr8UavzDykrK4tz587lu+/o0aP079+f6tWro1ar6dChA0FBQezcuROABg0a5OkdlJWVxalTp6hfv77Jfz8hiqpvkC+PP+qNjVqlb9npZlKBcrZqDl+O46mv/uTKLbmFaTbV2utuY9rYw6lNsHZI3v0yZZhFsshndFu3buX48eNkZ2eTmprKRx99RGRkJO3atcv3+JYtW/Lzzz9z+fJlFEXhwIED/PXXX/rncoMHD2bXrl1s27aNjIwMZs2ahZeXV4HlCWEOarWKBYMaM6d/EI0D3fFxc6BxoDtzBzRk71sdaVnNk9PXE+nz5X72nI01d7hlV9W2MGSTLtmdDoY1g3XbC1rVXJidRU4Btnz5cj7++GOuXbtGuXLlqF+/PtOnT6djx44A7Nu3jx49euhH52dlZTFlyhR++ukn4uPjqVSpEpMmTWLSpEn6Mjdu3Mibb75JZGQkjRs35ttvv6V27doGxyRTgAlzy87R8umvZ1my9yIqFbzS5RHGt6/O1hPXZYC5OVz+E75/AnIywa1ywT03SxFrrecsMtFZImu9AETpE3L8Om/8fIzUzBwqujpwKzkTrXLXpNFqFY8/6s2CQY0l2RW3y6GwvKduqjBXX3jttLkjeijWWs9Z5K1LIUTBejWoxOYJrang4kBsUgY5/yU5kAHmJe7mv7r/q9SQFAU/DgArTBSlnSQ6IUqhmt6uVHYvV+B+GWBeAu5+Jjdqh6435rlf4cf+kuwsjCQ6IUqp6wnpBe6TAebF7N6OJ35NYcQ2XbI7vwN+eEqSnQWRRCdEKVXYAHMVMsC8WOXX8aRyYxj1K9g5wYVdurF2kuwsgiQ6IUqpwgaYK4CTvY0MMC8uBU0ZVikIRv0GTp5wcBH88qYkOwsgiU6IUqqwAeZ2Nir2nbtJ7wX7OB4Zb8YoyyCf+jBsKzh56VYrD3kVtGVvAVdLIsMLDGSt3W5F6abVKgQfi2Llwctci0ulsrsTQ1oE0qKaB1M2nGDP2RvYqlVM7lyT8R2qY2sj321LTOwZWNEHUmKh8VDo/TkUMPOTpbDWek4SnYGs9QIQ1ktRFH48dIVZIadJy8qhob+GOQOCOBGZIAPMS8rNc7C8NyRHg39zGPELqG3yHnNkhe6ZX4cp5onxLtZaz0miM5C1XgDC+l26mcIra8MJvxqPWvXfoq4KMsC8pNy6AF93hIwEqNxU9wwvN9lZ2LRh1lrPWXY7Wgjx0Kp6OfPz2Jb0qOeDVtH1jZAB5iXIszq8uAcc3eHaYVjaBXKyLS7JWTNJdEKUAbY2amITCx53JwPMi5lHNXjhD3D0gKgw+KyGJLkSJIlOiDJCVjA3M/dAeHGvbtWDtDhwcIOgZ80dVZkgiU6IMuJBK5hn5miJKaTVJ0zgwi7QZoONA6TH657dZck5L26S6IQoIx60gvntlCw6fbaHb/ZeJCtHxn2Z3N3P5F4+Di7eEHsSvm4PWQW3tsXDk16XBrLW3kii7NBqFSatDmP7yRi02rzL+nSuU5EqHk58dyCCbK1CjYoufND3UVpU89SN05PhCA8nv44nyTdgSTvdqgeeNeHFP8De2ZxRWm09J4nOQNZ6AYiypaAB5rmJ63xsEu8Fn+TA+VsA+JQvx43EdF1vTWQ4QpHtmQ2ule7veJJyC5a0hcRrENAKnv8JHFzNEyPWW89JojOQtV4AQtxLURS2nYhm6sYTJKRl5XuMjVrFnP5B9GtUuYSjs0JpcbDyKV1vTL/HYPDPUK68WUKx1npOntEJIfJQqVT0alCJql4F30aT4Qgm5OgOQzfpklzkX/B9P13yEyYjiU4Ika/rCYUPR4iU4QimU648DNmgu30ZFfbfHJm3zB2V1ZBEJ4TI14OGI8SlZPHjoctkZksPTZNwcNXdtqzaDqJPwOI2ug4r9zqyQvfMTxhMntEZyFrvXQtRkE1Hr/HaumPkaO+/7lXcmUbMt3w5xnesQf+mftip1dJL82FlpcHXHeDGGd0QhBf3gquPbl8h04bpOxo9xLm31npOEp2BrPUCEKIghQ1HePxRb8a0q8aXu86z83QsAJXcHPB0ceB0dFK+x0svTSNkpcM3nXTj7Jwq6IYenN9ZaJIr7O/K0HNvrfWcJDoDWesFIERhHjQcAeBEZAKf/36OnadjCixHemkWQXYmLO0M0cd104ZpswucG7Ow1rcx595a6zlJdAay1gtACFN5fP5ezkQn5btPBTQOdGf9uFYlG1Rpl5MFn9XU9cK0sYfxB3WrIdzjqUUHCLsSn28Rxpx7a63nbM0dgBDCOsSlZha4TwHOxSSRkJZFeUc7/XZTPFeyVlqtwrEtiwhKSyABV9xzksj8qj22L/yOumItAGKT0tl09BrHIhMKLEcm7JYWncGs9ZuOEKby9H+tisL+lTjYqulVvxIDm/nTNNCdyWuOPvRzJWuk1SqsXvwhz8bM5Z3s0fyU0541dh/ymM1Z0lSO/NVpDSsuOPPHvzfyvV15N2nRSaIzmLVeAEKYSmHPidQqeKyqByciE0jJzAGggosDN1MyyO+fVVl/pnd00+c0OPoeU7NHszanIwAqtPxoN4tWNqdJURwYkPkeMc6P0K9hZdyd7Zm74195RlcASXQGstYLQAhTMaTnX1pWDiHHr7Pm7ysFPlMCeaa35tPxhMc7sea/JHeHwvd2s2lnc4IsOzdUQzZgG9BMel0+gCQ6A1nrBSCEKRnSSzNXkw93cCul4Od67k52hL7dmXJ2NnnLLqXP8wyJPyYxnX3nbvLuxhOkFzgQX2GG03qGazeA/X+DzANaGHXuC2LKei4kJIT//e9/nDhxAjs7O9q1a8f8+fPx8/MzSfnGkERnIEl0QpiWoc/0mlfzpG0NT/adu8mBC7dK5fO8AltcKhWNAzTUq1yeAxdu8m9M8gPL0rd26+6H3bPAzhmeW6ObUeUhmbKeW7VqFeXLl6d9+/aoVComTZrEmTNn+PPPP01SvjEk0RlIEp0QpvWgZ3r1fMtz6VYKSenZhZZT0DMoY1uAxXl8Yb/r3ap5OdOmphd2ahXL/7xMTj51Tp7fd/982PkeqG2h2Rjokc/UYEdWQNJ16DCl0M+G4q3njh8/TqNGjcjIyMDWtmQ7/MvwAiGEWfQN8mXHqehCnytpFYXwq/FMWnWU64np+ZaTo1V4f8tJIm6lULOiKzW9XQhwd+K1deF5yo5NzCA88hg7TkXf1wLMr8VV1ON/OxXN9N51uXgjhXOxyZyPTWbj0WuFJrlADyd+HNMcP3cnffnRiekFnpu+Qb66N7Z5GWzLwfa34NBXoORAz0/vFHz3lGFm9scff1CnTp0ST3IgLTqDSYtOCNMz9LlS8492EpOYYXC5KhX59uYEXWtxas86PNmoMi7lbHGwtTFoZpE+Qb4kZ2STkpHN5vBrfPrrWR7QQDOYj5sDB6d2ybPNqGduh7+Dra/o/tx0JPSeV+i8mAVRqfKW+9577zFjxoz7juvduzchISEFlnPp0iWqVKmif3306FE6duzIunXr6Nq1q0GxmJIkOgNJohPCfAp7nqcCHvF2ZWirQM7FJHMuNolDF2+TbWAWsrdRo1WUQo9XqzAqqTna2dDtUW9qVnShRkVXvvj9HKevJxYYv0l6mB79ETaP1/3ZvRrERxiV5MDwei4xMZHMzII7Enl4eKBW6xbHOXHiBF27duWzzz5j8ODBBsdiSnLrUghh8Ya0rEJ4ZAHP89QqxnWonucZ3YNagPY2Kmr5uJGckU1SejY3kwtvLSoK+Lk74uJgi4uDLccjE8jMKXh5ovKOtnz+bCP96/SsnIKfR6pVDGkRWOjnG6TR87ppwjaMhriLULmpUUnOGG5ubgYdd+LECbp06cLs2bPNluRAEp0QohR40PM8/TOr//hpHIlNzCiwBVXfT5OnBfWgFuO9La4HHV/5v2dtRY2/yLJSQaWGSg11C7geWVFsye5BTp48SZcuXZg5cyYjRowwSwy5ZOFVIYTFU6tVLBjUmDn9g2gc6I6PmwONA92Z0z8o36EFQ1pWKXC4QX4tqOI+3tj4i+TuZ3Iv7Nb9f+vLuu1m8Nlnn3Hjxg1eeeUVXFxc9D9Xrlwp8VjkGZ2B5BmdEKWHsTOFFPfxxa6gjidGdkix1npOEp2BrPUCEMJaGTtTSHEfX6z2zAbXSvknMwsZR2dOkugMZK0XgBBC5LLWek6e0QkhhLBqkuiEEEJYNUl0QgghrJokOiGEEFZNEp0QQgirJjOjGOHeCU+FEEJYPkl0BipNXW6ttYuwseQ8yDkAOQe5yvJ5kFuXQgghrJokOiGEEFZNEp0Veu+998wdgkWQ8yDnAOQc5CrL50GmABNCCGHVpEUnhBDCqkmiE0IIYdUk0QkhhLBqkuiEEEJYNUl0pdikSZPw9/fHzc2NypUr8/LLL5OZmQlAYmIizz33HG5ubnh7e/Phhx+aOdriU9B5iI2N5fnnn8fPzw83NzcaNWpEcHCwucMtFoVdC7liYmLw8PCgYcOG5gmymD3oHCxdupRatWrh7OxMlSpV2Lx5sxmjLT6FnYdTp07RuXNn3N3d8fHx4YUXXiA1NdXMEZcARZRap06dUpKTkxVFUZQbN24oHTp0UD788ENFURRl6NChSvfu3ZW4uDjl7Nmzir+/v7JixQpzhltsCjoPFy5cUD799FPl6tWrSk5OjhIcHKw4OTkpJ0+eNHPEplfYtZDrmWeeUTp16qQEBQWZIcLiV9g5WLJkiVK7dm0lLCxM0Wq1SnR0tHLhwgVzhltsCjsPQUFByrhx45SMjAwlNjZWeeyxx5QpU6aYM9wSIS26UqxOnTo4OzsDuinK1Go1586dIzU1lTVr1jBz5kw0Gg2PPPIIkyZN4ttvvzVzxMWjoPNQrVo1Xn/9dfz8/FCr1fTp04datWpx8OBBM0dsegWdg1ybN2/m9u3bDBkyxFwhFruCzkFOTg7Tp0/n888/p1GjRqhUKry9valWrZqZIy4ehV0LFy9eZPDgwdjb21OhQgX69u3LiRMnzBluiZBEV8rNnj0bFxcXKlasyLFjx5g0aRJnz54lMzMzzy2qhg0bcvz4cfMFWszyOw/3io2N5fTp0zRo0MAMERa/gs5BQkICr776KosXLzZzhMWvoH8PMTExhIWFUaVKFfz8/BgzZgyJiYnmDrfYFHQtvP7663z//fekpaURHR3Nxo0b6dOnj5mjLX6S6Eq5KVOmkJyczKlTpxg7diw+Pj4kJyfj7OyMre2dObs1Gg1JSUlmjLR45Xce7paZmcmzzz7LgAEDaNq0qZmiLF4FnYM333yT4cOHU7NmTTNHWPzyOwe3b98GYOfOnRw+fJjw8HAuXbrEK6+8YuZoi09B10KPHj3Yv38/rq6uVKpUCX9/f0aOHGnmaIufJDorUadOHYKCghg+fDguLi6kpqaSnZ2t35+QkICrq6sZIywZd5+HXJmZmTzzzDM4OTnxzTffmC+4EnL3Odi3bx8HDhzgrbfeMndYJerefw8Ab7/9Nl5eXnh5efH222+zZcsWM0dZ/O4+D3FxcXTp0oUxY8aQmprK7du3cXZ2ZvDgweYOs9jJMj1WJCsri3PnzlGrVi3s7Ow4duwYTZo0ASA8PJz69eubOcKSkXseQJfk+vfvT2ZmJps3b8be3t7M0ZWM3HPw+++/c/HiRXx9fQHIyMggLS0NLy8vTpw4QaVKlcwcafG5+99DuXLlzB2O2eSehwsXLpCWlsbkyZNRqVTY29vz4osv0qNHD3OHWPzM2hVGFFlSUpLy3XffKXFxcYpWq1WOHz+u1KlTRxkzZoyiKIoyZMgQpUePHkp8fLzy77//KgEBAVbZ67Kw85CZmak88cQTSufOnZW0tDRzh1psCjsHCQkJytWrV/U/c+fOVerWratcvXpVyc7ONnfoJvOgfw+jR49Wunbtqty+fVuJi4tTunbtqowePdrMUZteYechKSlJcXd3V7788kslKytLSUxMVIYMGaK0adPG3GEXO0l0pVRycrLSpUsXxcPDQ3F2dlaqVq2qvP7660pKSoqiKIqSkJCgPPvss4qLi4tSoUIF5f333zdzxMWjsPOwZ88eBVDKlSunODs7639mzZpl7rBN6kHXwt2WLVtmlcMLHnQOkpOTlWHDhinly5dXKlasqIwePVpJTEw0c9Sm96DzsH//fqV169ZK+fLlFQ8PD6VPnz5WO8zibrJ6gRBCCKsmnVGEEEJYNUl0QgghrJokOiGEEFZNEp0QQgirJolOCCGEVZNEJ4QQwqpJohNCCGHVJNEJIYSwapLohBBCWDVJdEJYuIyMDPr06YOPj0+e1RdiYmJo1aoV7du3p02bNvzzzz9mjFIIyyVTgAlh4TZt2sQff/zBjBkzaNu2rX4B3ZycHFQqFWq1mj179rB06VJ++OEHM0crhOWRFp0o06pUqcKmTZvMHYbemjVrGDBgQJ5tjzzyCGq1mszMzDxLLdnY2KBW6/4Jx8fHExQUVKKxGqtbt27s3LnT3GGIMkgSnSi1evbsycSJE+/bnpiYiJOTE7t27TJDVEWn1WqZOnUq06ZNy7O9bt26/PXXX1SsWJHOnTvn2Xfq1ClatWrFpEmTaN++fZ59GRkZaDQabt68+cDPLomE/8477/DGG28U62cIkR9JdKLUGjVqFKtWrSIjIyPP9tWrV1OpUiU6duxopsiKZtu2bXh4eNy3QO6VK1c4cOAAHh4erF69Os++unXr8ueff7J161YmTZqUZ9/u3bupX78+Xl5exR67Idq1a0d8fDwHDhwwdyiijJFEJ0qtvn37Ymtre19LZNmyZYwcORKVSkVMTAwDBgygQoUKBAQE8M4775CdnZ1veSqVivDwcP3r+fPn06FDB/3rKlWq8PHHH9OsWTOcnZ3p0aMHt2/fZvz48Wg0GmrWrMmff/6pPz45OZmJEycSEBBAxYoVGTp0KAkJCQX+PsHBwXTq1Om+7StXrqR8+fLMnTuXXbt2ERkZCZAnwWs0GpycnPK8b8uWLfTt21f/eu7cudSsWRNXV1eqV6/Ol19+CUD//v25cuUKgwYNwsXFhbFjxwI88NwZez5UKhWdOnUiODi4wHMgRHGQRCdKLTs7O4YMGcJ3332n33bq1CkOHz7M8OHDAXjuueews7Pj0qVL7Nu3j02bNvHJJ58U+TPXrl3Lhg0biIqK4urVq7Ro0YIuXbpw69YtnnvuOX2SABg5ciS3b9/m+PHjXLp0iaysrHxvteYKDw+ndu3a921fuXIlAwcOZODAgbi6uuo7nISFhdG+fXs6duzIsGHDmDNnTp733ZvoAgMD2bVrF4mJiSxdupQ33niDAwcOsG7dOgICAli9ejXJycksXrzY4HNnzPkAXQv07i8TQpQI8677KsTDOXnypKJWq5UrV64oiqIor732mtKzZ09FURQlMjJSAZTo6Gj98T/++KNSs2ZN/evAwEBl48aNiqIoCqAcPXpUv2/evHlK+/bt8xy7ePFi/es33nhDadGiRZ5YVCqVkpGRocTGxipqtVq5ffu2fv+///6r2NnZKdnZ2fn+LjVq1FDWrVuXZ9vBgwcVQDlw4ICiKIoyevRopU6dOg88L0ePHs3ze+bniSeeUGbOnKn/3XLPg6IYfu4MPR+5vv76a6VZs2YPjF8IU5IWnSjV6taty2OPPcaKFSvIzs7mhx9+YNSoUQBERkZSrlw5vL299cdXq1ZNf+uvKO4uy8nJ6b7XiqKQmppKREQEWq2WqlWrotFo0Gg0NGvWDLVaTXR0dL5lu7u7k5iYmGfbihUrqFmzJq1atQJg2LBhnD59mr///rvQOO9tzQH8+OOPNG7cGA8PDzQaDdu2bSuwo4qh587Q85ErMTERd3f3QmMXwtRszR2AEA9r1KhRzJ49m3r16qHVaunTpw8Afn5+pKenExMTo6+AIyIi8PPzy7ccZ2fnPJXy9evXixyTv78/arWaqKio+56dFaRhw4acOXNG/zozM5O1a9eSmpqKj49PnmO///57mjVrVmBZwcHBeW5lXrlyhWHDhrF9+3Y6dOiAra0t/fr1Q/lvGG3uMIVcxp47Q506dYqGDRs+VBlCGEtadKLUGzhwINHR0bzyyisMHToUOzs7ACpXrkzHjh15/fXXSUlJ4cqVK8yaNYthw4blW07jxo1ZuXIl2dnZhIeHs3LlyiLH5OPjQ79+/Zg4caK+1RQdHc3GjRsLfE+fPn3YvXu3/vWWLVtITk7myJEjhIeH638+/vhjVq9eTVZWVr7lXL9+nUuXLtG6dWv9tuTkZBRFoWLFiqjVarZt28Zvv/2m3+/t7c2FCxf0r409d4bavXs3vXv3fqgyhDCWJDpR6rm6ujJgwAAiIiL0ty1zrVq1irS0NAIDA2ndujW9evXizTffzLecBQsWEBoaikaj4a233nroSn358uX6W5Zubm60bduWI0eOFHh8z549uXnzpn4qr++//54BAwZQt25dfHx89D/jx48nMzOTkJCQfMvZunUrPXr0wMbGRr+tbt26vPPOO3Tq1AlPT0/Wrl2b59bm1KlT+fLLL9FoNIwfPx4w7twZYt++ffrzIERJkinAhLAgq1evZtOmTaxdu7bIZfTp04ehQ4fSv39/E0b28Lp3787rr79O165dzR2KKGMk0QlhZT755BPGjRuHq6uruUMRwiJIohNCCGHV5BmdEEIIqyaJTgghhFWTRCeEEMKqSaITQghh1STRCSGEsGqS6IQQQlg1SXRCCCGsmiQ6IYQQVk0SnRBCCKv2f5NvdRR/E8I/AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 405x315 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(4.5, 3.5), dpi=90)\n",
    "\n",
    "ax.plot(df.volume, df.energy, 'o-', label='energy')\n",
    "ax2 = ax.twinx()\n",
    "ax2.plot(df.volume, df.pressure, 'x-', color='C1', label='pressure')\n",
    "\n",
    "ax.set_xlabel('Volume (Å$^3$/atom)')\n",
    "ax.set_ylabel('Energy (eV/atom)')\n",
    "ax2.set_ylabel('Pressure (GPa)')\n",
    "\n",
    "fig.legend(loc='upper right', bbox_to_anchor=(0.85, 0.85));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f0ef8e3-f6d8-4d30-865b-e655d4cc5d5b",
   "metadata": {},
   "source": [
    "## Relaxation with a fixed cell\n",
    "\n",
    "In some scenarios, it is helpful to keep the cell constant but still relax the atomic configuration. This can be done by setting the `constant_cell` keyword argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "effd0748-0129-45ea-a1d0-4164d3b1d2d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "structure = bulk('PbTe', crystalstructure='rocksalt', a=6.6) \n",
    "structure.cell *= 1.1 ** (1 / 3)  # we strain the cell like above\n",
    "calc = CPUNEP('nep-PbTe.txt')\n",
    "structure.calc = calc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7900995-a5e5-4470-878e-640ab81f480b",
   "metadata": {},
   "source": [
    "We save the cell and atomic positions before relaxing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5c733775-db21-4007-852d-6e4a6875464a",
   "metadata": {},
   "outputs": [],
   "source": [
    "positions_pre_relax = structure.get_positions()\n",
    "cell_pre_relax = structure.get_cell()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6310a88-1a90-43fc-ae41-e38021b9522c",
   "metadata": {},
   "source": [
    "The structure is then relaxed under the constraint that the cell is unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ba77c2fc-be2c-4d05-a18a-ab0911da5fcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "relax_structure(structure, constant_cell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd868d9-fead-4a1f-aae5-a4265f38c148",
   "metadata": {},
   "source": [
    "Finally, we make sure that the atoms moved but that the cell did not change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d7483986-4955-431d-b1f8-aa976c2e54cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure that the atoms moved\n",
    "assert not np.allclose(positions_pre_relax, structure.get_positions())\n",
    "\n",
    "# And make sure that the cell didn't change\n",
    "assert np.allclose(cell_pre_relax, structure.get_cell())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
